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ABSTRACT 



A queuing system is investigated in which two types 
of customers. Type 1 and Type 2, attempt to access a multi- 
server operator switchboard. Customer arrivals for each 
type are assumed to be in accordance with independent 
Poisson processes with time constant parameters. Service 
times are assumed to be exponentially distributed and each 
type of customer has a different service rate. When a 
queue is formed. Type 1 customers are allowed to queue 

regardless of the length of the queue but Type 2 customers 
are blocked when the queue size exceeds a predesignated 
limit . 

An exact solution. of the problem was not obtained and 
a discussion of the complications which precluded an exact 
solution is presented. A simulation algorithm is developed 
and validated by comparison with an exact solution of the 
special case of the problem in which the service rates 
are equal . 

Finally, an analytical approximation is developed and 
shown to give reasonable results by comparison with the 
simulation model . 
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INTRODUCTION 



I . 

A. STATEMENT OF THE PROBLEM 

The problem addressed in this thesis arises in tele- 
communications systems and concerns the flow of, and delays 
to, customer traffic through a multiple operator-serviced 
switchboard. There are two types of customers. Type 1 and 
Type 2, and each type has its own arrival and service parameters. 
The arrival and service processes of the two types are sta- 
tistically independent. All switchboard operators can ser- 
vice either type of customer, and the queuing discipline is 
first-in, first-out (FIFO) . What makes this problem inter- 
esting is that while Type 1 customers are allowed to queue 
regardless of the queue size. Type 2 customers are blocked 
from joining the queue when the queue length equals or exceeds 
a predesignated limit. 

An example might be useful to ensure that the problem is 
understood. Consider the civil-sector telephone system, and 
let the Type 1 customers be those who are direct-dialing a 
collect, long-distance call, and the Type 2 customers be 
those who wish the operator to dial the call for them. The 
two types can be differentiated by the system on the basis 
of the number of digits they dial; the Type 1 customers would 
dial a total of eleven digits while the Type 2 customers 
would dial the single digit 0. Since the Type 1 customers 
require less work on the part of the operator, it may be 
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desirable under conditions of heavy use when the queue length 
gets large, to discourage customers from having the operator 
dial their calls for them by blocking Type 2 customers from 
joining the queue. In this case, both types of customers 
would be allowed to queue when the queue length is small, 
perhaps receiving a recording advising them that the operators 
are busy and asking them to wait. However, when the queue 
length exceeds its predesignated limit. Type 2 customers 
receive a busy sign'al or possibly a recording advising them 
to direct-dial their call. 

1. Simplifications and Assumptions 

In order to make the problem tractable, a number of 
complications which would exist in a real-world system have 
been ignored and a number of simplifying assumptions have 
been made. Among those elements of the problem which are 
not being considered are the following: 

a. Preemption by a high priority class of customers. 

A particular example would be customers who dial a three- 
digit code in the event of an emergency and who must receive 
preferential service. 

b. Balking or reneging. For purposes of this 
study, it is assumed that once a customer has joined the 
queue, he remains there until his service commences. 

c. The existence of a retrial population. In an 

actual situation, a Type 2 customer who is blocked from queuing 
would be expected to do one of three things: (i) decide not 

to attempt to complete the call at this time; (ii) try accessing 
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the system as a Type 1 customer; or (iii) continue trying 
to access the system as a Type 2 customer at repeated, random 
intervals. The repeated attempts to access the system by 
those Type 2 customers who exercise option (iii) can be seen 
to have an effect on the arrival rate, and these customers are 
said to belong to the retrial population. Although there 
exist a number of ways to model this phenomena, no attempt 
has been made to do so in this study. In addition, the 
increase in the arrival rate of Type 1 customers caused by 
those Type 2 customers who choose option (ii) is ignored. 

The assumptions which were considered necessary were: 

d. Arrivals are in accordance with a Poisson process 
and the arrival rates are constant over time for each type 

of customer. This is, admittedly, an unrealistic assumption 
since it is unlikely that the customer arrival rate will be 
the same at two A.M. as it is at two P.M. Although arrivals 
according to an inhomogeneous Poisson process can be modeled 
fairly well by computer simulation, the introduction of time- 
dependent arrival rates greatly decreases the chances of 
finding an analytical solution or approximation and it is 
for this reason that they have not been considered. 

e. Service times are exponentially distributed and 
service rates are constant. This assumption is somewhat less 
unrealistic since it seems plausible that service times should 
be the same, on the average, regardless of the time of day. 
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2 . Implications of Unequal Service Rates and the 

Queuing" Discipline 

The facts that the two types of customers have 
different service rates and that Type 2 customers are blocked 
from the queue when it becomes too long can be seen to have 
several important implications for the way in which the system 
functions. In particular, consider the following: 

a. Since the service rates are unequal, the system 
service rate is dependent not only on the number of custo- 
mers in service, but also on their types. Therefore, any 
solution technique, either analytical or simulation, must be 
able to keep track of the composition, by type, of the group 
of customers undergoing service at every time. 

b. When a change occurs in the state of the system, 
the new composition of the group of customers undergoing 
service depends not only on whether there are customers 
waiting in the queue, but also on the precise type of the 
customer at the head of the queue. Therefore, it is necessary 
to keep track of the composition, by type, of the queue, in 
detail. 

c. Since the blocking of Type 2 customers effectively 
sets their arrival rate to zero when the condition occurs, 

it can be seen that the overall system arrival rate is not 
constant, but is dependent upon the number in the system. 

This effect produces an interesting problem in computing 
the expected waiting time for an arriving customer. This 
problem is discussed in Chapter Three. 
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3 . Figures of Merit 



Three figures of merit were chosen to serve as indi- 
cators of system performance. They are: (1) the expected 

waiting time for a customer of either type who joins the 
system, (2) the proportion of time that there are no customers 
of either type in the system and (3) the proportion of time 
that the queue length equals or exceeds its limit and Type 2 
customers are blocked. For all three of the above, only the 
steady state or long-run value is considered. 

B. ANALYTICAL PROCEDURE 

1 . The Birth-Death Process 

The first attempt at a solution was to model the 
system as a multi-dimensional, birth-death process and then 
attempt to solve the resulting balance equations algebraically. 
The necessity of keeping track of the composition, by type, 
of both the group of customers undergoing service and the 
group of customers in the queue meant that the minimum number 
of dimensions which could be used to adequately model the 
state of the system was four. That is, the state variable which 
represented the total number of customers in the system would 
have four dimensions; one each for, the number of Type 1 
customers undergoing service, the number of Type 2 customers 
undergoing service, the number of Type 1 customers in the 
queue and the number of Type 2 customers in the queue. How- 
ever, even using a four dimensional state variable proved to 
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be inadequate to differentiate between the system state in 
which a Type 1 customer was first in the queue and the state 
in which a Type 2 customer was the first in the queue. Since 
the way in which the system changes states is dependent upon 
the type of the customer at the head of the queue , it would 
have been necessary to expand the state variable beyond four 
dimensions to model this effect. The complexity of the 
balance equations in the four dimensional case offered little 
hope of finding a solution and an expansion beyond four dimen- 
sions would reduce the likelihood of solving the equations 
still further. For this reason, the multi-dimensional, birth- 
death process approach was abandoned and other alternatives 
were considered. 

2 . Simulation 

The next attempt to model the system was via computer 
simulation since it was hoped that the development and exe- 
cution of the simulation model would provide insights into 
the workings of the system which would lead to ideas for 
analytical solutions or approximations. The development of 
the simulation model is covered in detail in Chapter two. 

3 . Verification of the Simulation Model 

Although verification of a simulation model is a 
normal step in the development of the model, it is mentioned 
separately here for several reasons. First, it was a major 
effort in itself; secondly, it utilized a generalization of 
Little's equation for the waiting time that might be of some 
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interest; and finally, it served as the basis for the first 
attempt at an analytical approximation. The verification 
of the simulation is covered in detail in Chapter Three. 

4 . An Analytical Approximation 

Chapter Four details the derivation of an approximate 
analytical solution based on the equations developed in 
Chapter Three for verifying the simulation model. In addition, 
the rationale for developing an analytical approximation 
model and the advantages to be gained from using such a model 
are presented. Finally, comparisons are made between results 
obtained from the analytical approximation and the simulation 
to determine the degree of agreement between the two models . 
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II. THE SIMULATION ALGORITHM 



A. RATIONALE FOR THE SIMULATION APPROACH 

Although it is generally accepted that an analytical 
solution is preferable to a simulation, there are a number 
of situations in which a simulation approach is indicated. 
Foremost among these is the case in which an analytical 
solution does not exist or, equivalently, an analytical solu- 
tion exists in theory, but is unobtainable in practice. As 
was mentioned in Chapter One, the latter situation is the case 
in the present problem. If it is impractical to obtain a 
complete and useable analytical solution, a search for a simple 
analytical approximation is likely to be worthwhile. Given 
the complexities of the problem of this thesis, derivation of 
an analytical approximation becomes a non-trivial task. 
Therefore, the simulation approach is utilized to gain insight 
into the problem and the operation of the queuing system in 
order to facilitate later analytical work. 

Another reason for simulating the system is that if an 
analytical approximation can be found, then the simulation 
can serve as a basis of comparison. This, of course, assumes 
a high degree of confidence in the simulation model which 
means that the simulation must somehow be verified. This 
important step in the development of the simulation model is 
covered in Chapter Three. 



14 



Finally, if a suitable analytical approximation cannot 
be found, then the simulation remains the only tool with 
which to study the system. 

B. CONCEPTUAL OPERATION OF THE SYSTEM 

The first step in the development of the simulation 
algorithm was to study the operation of the system so as 
to make sure that the manner in which the system changed 
from one state to another was fully understood. The system 
was viewed as staying in a current state for a random time 
after which one of four statistically independent events would 
occur with certain specified probabilities. The occurrence 
of these events causes the system to change state. The four 
possible events are: (i) a Type 1 customer arrives and joins 
the system, (ii) a Type 2 customer arrives and may or may 
not join the system depending on the current state, (iii) a 
Type 1 customer completes service, or (iv) a Type 2 customer 
completes service. The conceptual operation of the system 
and the four events along with their probabilities of 
occurrence are displayed in Figure 2.1. Figure 2.2 shows the 
system changes associated with the occurrence of each of the 
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T = queue threshold (applies to Type 2 customers) . 

N = total number of customers in the system 

= J1 + J2 + Q1 + Q2 

A^ = Type 1 arrival rate. 

A 2 = Type 2 arrival rate. 

y^ = Type 1 service rate. 

= Type 2 service rate. 

d = A^ + A 2 + Jly 1 + J2y2 

A few notes will be helpful to fully understand the two 

figures . 



1. The "sojourn time" is the random time period 
during which the system stays in its current state. The 
sojourn time is exponentially distributed with mean equal to 



A^ + A 2 + Jly^ + J2y2 



which is j. Note that this means that 



the distribution of the sojourn time is not fixed, but that 
it changes with each change in J1 and/or J2. 

2. According to the queuing discipline assumed for 
this problem, a Type 2 customer can only join the system 
when N, the number in the system, is less than T, the queue 
threshold. This is illustrated in Figure 2.2, where it can 
be seen that no state change takes place if a Type 2 arrival 
occurs when the number in the system is larger than the cutoff 
value . 



3. Figure 2.2 also illustrates the fact that when 
a queue is formed, the new state of the system depends upon 
whether a Type 1 or Type 2 customer is at the head of the 
queue. 
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FIGURE 2.1. Possible Events and Their Probabilities of Occurrence 
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FIGURE 2.2 
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C. STRUCTURE OF THE ALGORITHM 



After the operation of the system was understood and 
represented by a probability model, the next step was to 
write out a simulation algorithm in flow-chart form. The 
algorithm must be capable of performing at least the follow- 
ing functions: (i) it must automatically choose an appro- 

priate value for each sojourn time based on its current proba- 
bility distribution, (ii) it must determine which of the 

/ 

four possible events is to occur next, (iii) it must update 
the state of the system based on the type of event that 
occurred, the value of Jl, J2, Ql, Q2, and the type of cus- 
tomer at the head of the queue, if formed, and finally 
(iv) it must gather appropriate statistics on the figures 
of merit. The general flow-chart of Figure 2.3 describes 
the way in which the program performs these four functions. 
Notice that the algorithm essentially consists of one loop 
embedded in another. The inner loop performs the major 
tasks of choosing events and sojourn times and updating the 
system variables. This loop continues to iterate for a 
specified length of simulation time, for example one hour. 

When the time limit is reached, control is passed to the outer 
DO loop which gathers interim statistics and then checks to 
see if an iteration counter has reached its maximum value. 

If it has not, then control is passed back to the inner loop 
and the process starts anew with the values of the system 
variables preserved from the end of the previous iteration. 



19 



GENERAL FLOW CHART FOR SIMULATION ALGORITHM 
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FIGURE 2.3 
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If the iteration counter has reached its upper limit, then 
final statistics are gathered and the algorithm terminates. 

In this manner one day of system operation can be simulated 
by setting the upper limit for the inner loop to one hour 
and the upper limit for the outer loop to twenty-four. The 
algorithm uses four sub-algorithms called R0UTE1 , ROUTE 2 , 

R0UTE3 and.R0UTE4. The sub-algorithms are utilized when an 
event occurs and their function is to correctly update the 
system variables according to the rules depicted in Figure 2.2. 
Therefore, for instance, R0UTE2, which is utilized when a 
Type 2 arrival occurs , would first check whether all the 
servers were busy; if not, it would increase the number of 
Type 2 customers in service by one. If the servers were all 
busy, R0UTE2 would then test whether the number in the system 
equaled or exceeded the queue threshold; if not then the 
number of Type 2 customers in the queue would be increased 
by one. If the number in the system were greater than the 
queue threshold, there would be no change in the state of 
the system and control would be passed back to the main 
algorithm. The flow-chart for R0UTE2 is shown in Figure 2.4. 
The flow-charts for the other sub-algorithms are similar. 

The composition of the queue is maintained in a vector 

th 

of one's and two's, called NEXT1, in which a one in the i 

Vi 

position indicates that the i customer in the queue is 
a Type 1. Each time a customer arrives, another element is 
added to the vector; either a one or a two depending on 
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FIGURE 2.4. 

General Flow Chart 
for ROUTE 2 
Subroutine 
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the type of customer who arrived. Each time a service event 
occurs, a counter is incremented by one. When this counter 
is used as the subscript of NEXTl, it points to the element 
the vector which represents the customer currently at 
the head of the queue. For example, assume that five arrival 
events have occurred and that the composition of the vector 
NEXTl is as shown in Figure 2.5. Let the value of the 
service counter be four indicating that the fourth service 
event is about to occur. Then NEXTl (4) indicates that the 
customer at the head of the queue is a Type 2. 

Once the problem was understood and the flow-chart 
written, the next step was to choose an appropriate computer 
language and program the algorithm. 

D. THE COMPUTER PROGRAM 

Fortran IV was chosen as the language to implement the 
simulation algorithm. Initially, consideration was given to 
simulation-oriented languages such as SIMSCRIPT and GPSS, 
but it was felt that their sophistication and abstraction 
would tend to obscure some of the inner workings of the 
queuing system and would thus offer less insight into the 
problem. Coding of the program was a straightforward trans- 
lation of the flow-chart. A program listing, with explana- 
tory comments and a sample output, is provided at the end 
of the thesis. 

The three figures of merit are printed out after each 
hour. These are cumulative statistics - that is, they 
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Vector NEXTl 
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'\ 




NEXTl (4) points to the fourth element in NEXTl 
and indicates that when the fourth service event 
occurs, the customer at the head of the queue is 
Type 2 . 



FIGURE 2.5. Method of Identifying the Type of the 

Customer Currently at the Head of the Queue 



25 



reflect the operation of the system from the start of the 
simulation, and not just from the previous hour. The printing 
out of cumulative statistics at the end of each hour was 
necessary in order to ensure that the simulation had been 
allowed to run long enough to reach the practical equivalent 
of steady— state . Note that there is little fluctuation in 
the values of the figures of merit toward the end of the 
program run. If there had continued to be significant fluc- 
tuation, it would have indicated that the system was not 
in steady-state and that the program should be run again for 
longer simulation time. The program keeps track of the 
composition of the system, by customer type, at all times 
and, for each customer who joins the queue, it computes the 
actual time he waits for service. Thus, the average waiting 
times are computed by simply totaling all of the waiting 
times and dividing by the number of customers who joined 
the queue. The probability that there are no customers in 
the system is actually the long-run fraction of time that 
the system was in that state. Similarly, the probability 
that a Type 2 customer is blocked is the long-run fraction 
of time that the queue threshold was exceeded. 

The program utilizes the IMSL sub-routine RANDOM to 

r 

generate the pseudo-random numbers needed for determining 
the sojourn time and choosing the next event. RANDOM takes 
as input a "seed" number and, for a given seed, will always 
produce the same string of numbers which are uniformly 
distributed on the range from zero to one. 
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Looking at the print-out from the simulation program, 
the first question that might be asked is: "Are the 

results meaningful — do they accurately represent the 
functioning of the real system?" Thus the need for some 
kind of verification of the simulation program becomes 
obvious. The process of verification is detailed in the 
next chapter. 
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HI- VERIFICATION OF THE SIMULATION ALGORITHM 



A. NEED FOR VERIFICATION 

Verification of a simulation model by testing it 
against the results obtained from a mathematical model is 
the final step in its development. As mentioned in Chapter 
If the verification process has been separated from the 
rest of the development of the algorithm in order to better 
highlight some interesting results. Of course, as is often 
the case, the reason for going to a simulation approach in 
the first place was that no general mathematical model 
existed. However, it was felt that if a mathematical model 
could be developed for a special case and the simulation 
verified against it, then it would at least indicate that 
the simulation was performing properly. This is the approach 
that was taken, and it is detailed below. 

B. METHOD OF VERIFICATION 

1 . The Mathematical Model 

It was decided to attempt to model the special case 
of the problem in which the service times of the Type 1 and 
Type 2 customers are equal. All other aspects of the prob- 
lem remain unchanged; that is, the arrivals occur according 
to a homogeneous Poisson process with different parameters 
for each type of customer, and Type 2 customers are blocked 
from joining the queue when the queue length equals or 
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exceeds a given threshold. Comparison between the mathe- 
matical model and the simulation model was made on the basis 
of the three figures of merit mentioned in Chapter I; those 
were, (i) the probability that the system was empty, (ii) the 
expected waiting time for any customer and (iii) the probab- 
ility that a Type 2 customer would be blocked. The equations 
for (i) and (iii) are, although complicated and tedious, 
straightforward solutions of the balance equations which 
resulted from application of the theory of birth-death 
processes.^" The derivation of these equations is detailed 
in sections two and four of this chapter. 

It appeared that the simplest way of computing the 

2 

expected waiting time was to utilize Little's equation ; 



E [W] 



E [Q] 

X 



(3.1) 



Where : 



E [W] = the expected waiting time of any 

customer 

E [Q] = the expected number of customers in 

the queue 

X = the system arrival rate 



^Gaver, D.P., and Thompson, G.L., Programming and 
Probability Models in Operations Research , p. 467-482, 
Brooks/Cole, 1973. 

2 Hillier, F.S., and Lieberman, G.J., Introduction to 
Operations Research, p. 291—292, Holden— Day , 1967. 
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Little's equation, as shown above, assumes a constant 
arrival rate^ and is not directly applicable to the problem 
under study since, as was discussed in Chapter I, the 
system arrival rate depends upon the number of customers 
in the system. Therefore, it was decided to modify 
Little's equation as follows: 



E[W] = 

X E 



(3.2) 



Where : 



Ag = the "effective" rate at which customers 
join the system. 



A 



E 




T-l 



l P 

n=0 



n 



(3.3) 



p = the probability that there are exactly 
n customers in the system. 

Before equation (3.2) can be verified, it will be 
necessary to set up the balance equations and derive an 
expression for p n . At the same time, an expression for the 
probability that the system is empty will be derived. 



3 



Hillier and Lieberman , 



P. 



291-292. 
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2 . A Derivation of an Expression for the Probability 
That the System is Empty 

The theory of birth-death processes provides 
expressions for long-run or steady-state probability that 
the number of customers in the system is a given value. 

To determine these probabilities, it is first necessary 
to define some terms. Let 



N 



P 

X 

X 



n 

1 

2 



u 



X 



(i) 



T 



s 



u 



(i) 



the number of customers in the system, 

the probability that N = n, 

the arrival rate for Type 1 customers, 

the arrival rate for Type 2 customers, 

the service rate for both types of customers, 

overall system arrival rate when N = i, 

queue threshold; that is when N >_ T, Type 2 
customers are blocked from the queue, 

the number of servers, 

system service rate when N = i. 



It can then be shown that an expression for 




, (0) .(1) , (2) . (n-1) 

X • X • X • • • X „ 

P n " (1) (2) (3) (n) p 0 
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Gaver and Thompson, pp. 



467-482. 
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which comes directly from the basic balance equation, 



This leads to the following expression for p 



0 



00 i , (j-1) 

1 + l { n } 

i-1 j=l y(D) 



(3.5; 



The remainder of this section is devoted to manipulating 
equation (3.4) into a form suitable for use with the problem 
at hand. 

The first step is to write down the system arrival 
and service rates that result from the queuing discipline: 



(0) = a + A 

A x + A 2 


(0) 

y 


(1) = A + A 

X 1 A 2 


(1) 

y 




(2) 



= 0 



= y 



= 2y 



X (T " 1} = A, + * 2 



X (T) = A. 



(s) _ 



= sy 



(s+1) _ 



= su 



A (n) = A. 



y (n) = sy . 
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Rewriting equation (3.5) in terms of the above expression 
yields an expression for ^ in terms of the system service rate 
and the rate at which customers actually j pin the system . 

This is an important point and it will be repeatedly 
emphasized in the remainder the the thesis. In a "normal" 
queuing system, the arrival rate is identical to the rate 
at which customers join the system. However, a distinction 
must be made in any system in which the queuing discipline 
allows customers to be blocked from the system. The correct 
expression for p^ , in terms of the rate at which customers 
join the system, is shown below: 



s 

1+ [ y 

n=l 



(x 1+ a 2 ) 



n 



n : u 



n 



T 

I 



Ul+x 2 ) 



n 



n _ c ,i (n-s) n 

n-s+l sis u 



, .T (n-T) 

y <> l +> 2 ) *1 

n=T+l s!s (n - s) u n 



1 



Next examine each of the summations in the 
denominator to see if they can be summed in closed form. 

s (A.+A ) n 

a. £ — . This expression cannot be summed 

n=l ni w 

in closed form; however, the sum can easily be computed in 
two ways; (i) by direct summation if s is not large or 
(ii) by reference to a table of cumulative Poisson proba- 
bilities. The second method is possible since 



( (A^+A 2 )/ u) 

nl 



n 



exp (- 



X l +X 2, 
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is the density function, with respect to counting measure, 

A,+A 9 

for a Poisson distribution with parameter — . Thus 

y 

if the figure obtained from a cumulative Poisson table 

look-up with the appropriate parameter is dived by 
X l +X 2 

exp(- — ) then the desired sum will result. 

T (A 1 + A 2 ) n 

h. 1 / ■ . — — . This sum can be expressed 

n=s+l s ! s u 



in closed form as follows: 



■i 00 (A,+A~ 

Tft l 12 



n 



n 



00 ( A, +A 0 ) 

1 ■ vk- T- . l 



s! ' L , , (n-s) n ^ , .. (n-s ) n 

n=s+l s y n=T+l s y 



s+1 

•i (A,+A 9 ) 00 (A, + A ~ ) ■ 

l r 1 Z v t 1 z \ 3 

Si 1 S + l Sy ; 

Sy 3 = 0 



( A , +A 2 ) (T+1) » A , +A 2 , 

s (T + l-s ) ^ (T+l) ] ] 



( W 



(s+i: 



( A ! + A 2 ) 



(T+l) 



s! (1- 



X l +X 2 , 

sy J 



sy 



(s+1) 



(T+l-s) (T+i: 

s y 



oo M ( n_T ) 

y (X 1 +X 2’ X 1 

n=T+l sls (n - s \ n 



This sum can also be 



expressed in closed form: 
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- (X 1 + X 2 ) T A 1 < n -T) 

n=T+l s!s ln_s) i, n 



( A A 2 ) 00 



(n-T) 



s ! 



1 



n=T+l s ln - s) „ n 



(\ 1 +X 2 ) x 1 

- — ( ± ) ( — - ) 

s! k o (T+l-s) (T+l ) M A , 
s U 1 1 

~Sy 



Combining the above sums yields the desired 



expression for : 



n 



s ( A , +A 0 ) , 

1+ [ 1 I + (A-<- 



n=l n ! 



n s ! 
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A i + A 2 



U l +X 2 ) 



(T+l) 
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The next step is to derive equation (3.2), 
the modified form of Little's equation. 

3 . Derivation of Little's Equation for the Case Where 
the Service Rate is a Function of the Number of 
Customers in the System 

The validation of equation (3.2) proceeded in two 
steps; first, the equation was derived in general to 
be applicable to any queing system and next, the general 
equation was modified to yield a form directly applicable 
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to the problem at hand. Before beginning the derivation, 
it was necessary to define a number of quantities . 



T = the queue threshold. 

A (n) = combined rate at which customers join the 

system when N = n. 

A = the event that a customer, of either type, 
arrives and joins the system. 



W = waiting time for any customer . 
P n = probability that N = n. 



In 

expression 



the general case, the desired form of the 
for the expected waiting time is : 



total expected waiting time for all customers 

_ who join the system during time t 

total expected number of customers who join 
the system during time t. 



In order to write the above expression in mathe- 
matical form, first consider the expected waiting time of 
a customer, given that a customer arrives and joins the 
system when N = n. This can be seen to equal zero if there 
are any idle servers and to equal the expected service 
times of all customers ahead of him in the queue plus the 
expected service time of the customer being served if the 
servers are all busy and a queue is formed. This can be 
expressed as: 
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0 



if n < s 



E [W [ A,N=n] = - 






n-s+1 

s 



if n > s 



Now consider the expected contribution to the total 
expected steady-state waiting time which results from a 
customer arrival in a short time interval (dt) . This is 
simply the expression derived above, multiplied by the 
steady-state probability that a customer joins the system 
during a time interval dt when there are n customers in the 
system (A(n)dt), and by the probability that there are n 
customers in the system (p n ) , and then summed over the 
range of N. The expression appears as: 



{ n - - S ~— A (n) p } dt . 

L su n 

n=s 

If this last expression is integrated over the 
range (0,t) then, in the limit as t goes to infinity, the 
result is equal to the total expected waiting time for all 
customers who joined the system during time t, which is the 
numerator of equation (3.7). That is: 



t 



/ 

0 



oo 



l i 

n=s 



n-s+1 

sy 



•A (n) P n >dt 



oo 



l 



. n-s+1 

1 Sy 



A(n)p n >t . 
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The denominator of equation (3.7) is simply: 



00 



l (A(n)p )t . 
n=0 n 



Combining the expression for the numerator and 
denominator yields the desired expression for the expected 
waiting time which is applicable in general: 



E [W] 



l ( 



n=s 



n-s+1 

sy 



■A (n) p n > 



l A (n) p 
n=0 n 



Equation (3.8) can be used to determine the expected 
waiting time at steady-state for any queuing system in 
which the rate at which customers join the system is a 
function of the number of customers in the system. In order 
to make equation (3.8) more directly applicable to the 
particular system under study, it is necessary to define 
some more terms. They are: 

a:, (n) = the probability that a Type 1 customer 

can join the system upon arrival when 
there are n customers in the system. 

id- (n) = the probability that a Type 2 customer 

can join the system upon arrival when 
there are n customers in the system. 
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Note again the important distinction made between 
the rate at which customers arrive ( x ^ and a ) and the 
rate at which customers joint the system by either join 
directly into service or joining the queue (A(n)). The 
two rates are related as follows: 

A(n) = A lUi (n) + A^tn) (3.9) 

Combining equations (3.8) and (3.9) yields the 
desired expression for the expected waiting time. 



E [W] 



i l ' 



n=s 



n-s+1 

sy 



’l (n) P n } + X 2 l {S if ±1 “2 (n) Pn 1 



n=s 



X 1 ^ “l (n,p n +X 2 I “2 (n)p n 
n=0 n=0 



(3.10) 



Equation (3.10) can now be applied directly to the 
problem of interest by inserting appropriate values for 
A^,A 2 >w^(n) and (n) . However, in order to display the 
connection between equation (3.10) and equation (3.2), the 
modified form of Little's equation, some additional work 
is necessary. 

Recall equation (3.4) in section two of this chapter: 

, (0) .(1) , (2) . (n-1) 

A • A • A • • • A ^ a 

p n " TT5 (21 (Tj (771 p o • u 

y • y • y • • • y 
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Writing equation (3.4) in recursive form yields; 



P 



n 



x (n) p n+l 



(3.11) 



Now apply the definitions above to substitute for A ^ ; 



s y p 



n+1 



n 



A^a)^ (n) + A^w^Tn) 



Now insert this expression for p^ into the numerator 



of equation (3.10); 



A, l ( ^ — 0 ). (n) p } + X- l ( u ~— -a), (n) p } 
1 L Su 1 n 2 L s y 2 n 
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S y P 
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n=s 



(n-s+l) (n) P n+1 

A 1 (u 1 (n) + A 2 oj 2 (n) 
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Letting m - n+1 yields the following expression 
for the numerator of equation (3.10): 



A 



1 



oo 



l 

m=s+l 



(m-s) id. (m-1) p 
± m 

A ^ oj ^ (m-1) + A 2^2 (m-1 ) 
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l 
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1 m 

A^w^ (m-1 ) + A 2 ^ 2 (m-1) 



oo 
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m=s+l 



A^o)^ (m-1 ) + A 2 U 2 ( m “l) 

Aw (m-1) + Xu (m-1) (m-s)p i 
11 " “ 



oo 

= l 

m=s+l 



(m-s)p m 



(3.12) 



Notice that equation (3.11) is an expression for 
the expected number of customers in the queue in steady- 
state. This allows us to rewrite equation (3.9) as follows: 



E [W] 



E[Q] 



1 ^ u l (n)P n 

n=0 



+ X I “2 (n)p n 
n=0 



expected number in the queue 

(expected value of the rate at which 
customers join the system) 



Thus, equation (3.10) can be seen as a generalization 
of Little's equation (3.2), as would be anticipated. To 
complete the derivation it is only necessary to insert values 
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for w^(n) and (n) as appropriate for the particular 
problem. In the present case, these values are 



(n) 



1.0 for all values of n 






1.0 if 0<n<T 




0.0 if n > T 



Inserting these values into equation (3.9) yields 



E [W] 



E[Q] 



T-l 




which is identical to equation (3.2). 

4 . Calculations of the Expected Waiting Time and the 
Probability That a Type 2 Customer is Blocked 

Having derived equation (3.2) as the expression for 

the expected waiting time, it now remains to obtain 

expressions for the expected number in the queue ( [EQ] ) and 

the expected rate at which customers actually join the 



system (X^ + ^ I P n ) • Once again, these calculations 

n=0 

are straightforward but tedious and, hence, they are presented 
in outline form. First consider the following expression 
for the expected number of customers in the queue: 



T-l 
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E[Q] 
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= l (n-s)p 
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n=s 



n 



Both of the summations above can be expressed in closed 
form. Details of the necessary calculations are included 
in Appendix A and only the results of those calculations 
given below: 
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Next, consider the expression for the expected rate 
at which customers actually join the system; 
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“ \ + X 2 *P„ 
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Both of the summations above were encountered in 
deriving the expression for p^; the first sum cannot be 
expressed in closed form, but the second can and details of 
the calculations can be found in section two of this chapter. 
The final form for A is: 
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The final figure of merit is the probability that 
a Type 2 customer is blocked from joining the queue (P,^) ' 
This is simply the probability that the number in the 
system equals or exceeds the queue threshold which is: 
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Since the summation, j p was just encountered in finding 

n=0 n 

Ag it is possible to write the desired expression for 
directly : 
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5 . Results of the Validation Efforts 

At this point, expressions have been derived for all 
three of the figures of merit; the probability that the system 
is empty (p^) , the expected waiting time for any customer 
(E[W]), and the probability that a Type 2 customer will be 
blocked from joining the system ( P 2 j_j ) • The next step is to 
compute values for these quantities and compare the results 
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obtained with the values of those same quantities which 
result from running the computer simulation. 

In order to reduce the variance in the simulation 
results, the method of antithetic variables was used. 5 
Recall from Chapter Two that the RANDOM subroutine utilized 
by the program produces a string of pseudo-random numbers 
which range from zero to one and that it will produce the 
same string each time unless the "seed" number is changed. 

In the antithetic variable method of variance reduction two 
simulation runs are made for each data point desired. In 
the first run the random numbers are used, unmodified, to 
determine the sojourn time and choose the next event. In the 
second run, the program is modified slightly so that the 
random numbers are subtracted from one forming an antithetic 
random variable which is then used as before. Results from 
the two runs are then averaged and the value becomes the 
data point for that iteration. The seed number is then 
changed to produce a new random number string and the process 
is repeated for as many iterations as desired. Finally, the 
data points for all of the iterations are averaged to produce 
an estimate of the variable under consideration. As illus- 
trated in Gaver and Thompson [Ref. 4], sampling by antithetic 
variables can reduce the variance in simulation results by 
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Gaver and Thompson, 



pp. 



584-586. 
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more than one-half of the variance expected using straight- 
forward sampling techniques. 

Finally, the comparison of simulation' and analytical 
values, which serves as the validation of the simulation model, 
is presented in Figure 3.1. The values shown under the 
simulation heading were obtained by averaging over the results 
of ten iterations of the antithetic variable sampling tech- 
nique just described; that is, twenty simulation runs were 
made in all. Each run simulated thirty-six hours of system 
operation. The comparison, measured in terms of the percent 
difference between the analytical and simulation values, 
seems to offer a good indication that the simulation model 
is functioning as intended. Notice that the largest discrepancy 
is consistently in the expected waiting time values. No 
explanation for this evidently systematic effect was immediately 
apparent. 
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Verification of the Simulation 





Case I A^ = 1.0 
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IV. AN ANALYTICAL APPROXIMATION 



A. RATIONALE FOR DERIVING AN ANALYTICAL APPROXIMATION 

In the previous three chapters a simulation model was 
developed and verified. It was seen that the simulation 
appeared to do a reasonably good job of modeling the system; 
good enough, at least, for most engineering purposes. Why, 
then, go to the trouble of developing a second model? There 
are several reasons. First, all of the reasons for developing 
the simulation apply equally to an analytical model. In 
particular, the process of developing an analytical approxi- 
mation allows an opportunity to gain additional insight into 
the problem from a different perspective. Such insight could 
lead further toward the ultimate goal, which is an exact 
but simple and useable analytical solution if it exists. 
Secondly, the implementation of the approximation described 
in this chapter requires tools no more sophisticated than a 
pocket calculator. Thus, the analytical approach offers 
some practical, financial advantages over the simulation. 
Financial advantages are even more apparent when the need to 
do sensitivity analysis is considered. In order to determine 
the sensitivity of the system to a change in one variable, 
it is necessary to do a large number of simulation runs for 
each different value of the variable. Consider Figure 3.4 
in which it was necessary to do eighty simulation runs in 
order to be able to compare four variable changes. The 
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analytical approximation need only be done once for each 
variable change and if done on a computer, the calculations 
required take only a fraction of the time used by one simula- 
tion run. On the other hand, it would be easier to modify 
the simulation to investigate other figures of merit than 
it would be to write new equations for the approximation 
model. 

B. DEVELOPMENT OF THE ANALYTICAL APPROXIMATION 

Recall that in Chapter Three it was shown that an 
analytical solution to the problem existed in the special 
case in which the service rates were equal and that expressions 
were developed for the expected waiting time, the probability 
that the system is empty, and the probability that a Type 2 
customer will be blocked. This suggested approximating the 
unequal service rate case by taking a weighted average of 
the service rates and using the resulting value in the equa- 
tions for the equal service rate case. The intended weighting 
was 
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Where : 



y = the weighted average of service rates. 



50 



All other variables are as defined previously. The weighting 
factors, which are the ratio of each type of customer's 
arrival rate to the overall system arrival rate, are intended 
to represent the long-run proportion of time that each 
operator would spend serving each type of customer. Unfor- 
tunately, this weighting scheme is infeasible since it 
requires values for the p n 's which are not available. There- 
fore, an iterative procedure was decided upon in which another 
type of weighting was used to get approximate values for the 
P n ' s which could then be used in equation 4.1 to get better 
approximations which are again used in equation 4.1 and so 
on. The initial weighting utilized was 



H' 



M 1 A. +X~ y 2 A, +A. 



(4.2) 



The rationale behind this weighting scheme is the same as 
before, however, it can be seen that this weighting does 
not take into account the effect of blocking Type 2 customers. 

A computer program was written which would perform the 
necessary calculations and iterate the procedure as many 
times as desired. A program listing and sample output are 
provided at the end of the thesis. 

C. VERIFICATION OF THE ANALYTICAL APPROXIMATION 

As before, the final step in developing a model is to 
test its performance by comparison with a suitable mathe- 
matical model. Recall that the simulation model was verified 
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by testing it against an analytical solution to the special 
case of the problem in which the two types of customers have 
equal service rates. This would suggest trying to find another 
special case of the problem in which the service rates are 
different and for which an exact solution could be found 
with which to verify the analytical approximation. Unfor- 
tunately, such a special case was not discovered. An 
attempt was made to solve the problem for the case of a 
single server with the queue threshold equal to one where 
Type 2 customers could only access the system when it was 
empty. Although it seemed that a solution should be obtain- 
able, the complications associated with the different service 
rates made even this comparatively simple problem unsolvable 
in the time available. 

Since the approximation model cannot be verified analytically, 
the only alternative available is to verify it by comparison 
to the simulation model. Of course, the simulation is itself 
an approximation since it is only exact in the limit as time 
goes to infinity and true steady-state is reached. However, 
it was shown in Chapter Three that, at least for small systems 
and equal service rates, the results obtained from running the 
program for thirty-six hours of simulation time differ from 
the exact analytical solution by less than three percent. 

Chapter Two described how the simulation program was able 
to keep track of the identity of the customer at the head 
of the queue and thereby correctly adjust the system service 
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rate when a state change occurred. Therefore, the simula- 
tion should be equally valid in the case where the service 
rates are different. Thus, a reasonable amount of agree- 
ment between the approximation and the simulation, while not 
conclusive, should be sufficient to provide confidence in 
the approximation model. 

To see how closely the two models agree in two test 
cases, refer to Figure 4.1. The same method of antithetic 
variable sampling that was discussed in Chapter Three was 
used to obtain the results for the simulation model. The 
approximation results were obtained by carrying out ten 
iterations of the procedures discussed earlier in this chap- 
ter. It is interesting to note that the largest change in 

the figures of merit occurred between the first iteration — 

A X 2 

which utilized y' = y^ ^ — + y 2 ^ and the second 

iteration, and that there was little change between succeeding 
iterations. In fact, there was no change at all after the 
fifth iteration. This effect is apparent in the sample 
output. 

The degree of agreement between the two models shown in 
Figure 4.1 indicates that either model would be appropriate 
for a first estimate of system performance. The choice of 
model would therefore depend on the relative merits of each 
as evaluated in a particular situation. For instance, if 
sensitivity analysis was to be performed, the approximation 
model would be appropriate, while if an estimate of the 
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Analytical 

Approximation 



Simulation 



Percent Difference 
( (Sim-Anal) / (Sim) 



Case I A^ = 1.0 A 2 

s = 

E[W] .45110463 
p Q .1047012 

P 2b -10048 



•5 y ^ = .666 


U 2 = .833 


T = 6 


. 4700145 


4.02 


. 099507 


5.22 


. 1064845 


5. 63 



Case II \ = 1.75 A 2 = .249 v 1 = 1.0 U 2 = .833 

s = 3 T = 6 



E [W] 

P 0 

p 2b 



. 39178724 
. 10854438 
.11799 



. 40591 
. 109263 
.12258 



3.47 
. 657 
3.74 



Figure 4.1 
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proportion of customers whose waiting time was greater than 
a given value was desired then the simulation model should 
be chosen. 
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V. CONCLUSION 



Two aspects of this problem which made it interesting 
were the different service rates for the two types of 
customers and the fact that Type 2 customers are blocked 
from the system whenever the queue threshold is reached. 
Recall that Chapters One and Two discussed the implications 
of these aspects of the problem and showed that they made 
an exact analytical solution difficult, if not impossible, 
to achieve. In particular it was shown that the different 
service rates required that any solution be able to keep 
track of both the number and type of customers in the queue 
and in service as well as the identity of the customer 
at the head of the queue. Also, the blocking of Type 2 
customers implies that the systems arrival rate is not con- 
stant but is a function of the number of customers in the 
system. 

The only real attempt made to find an analytical solu- 
tion utilized the theory of birth and death processes and 
it is possible that some other analytical approach could 
yield results. However, given the complications mentioned 
above, it seems that any future effort expended on this 
problem could best be directed toward improving the approxi- 
mate solutions. 

It would be worthwhile to attempt to verify both the 
simulation and the analytical approximation with an exact 
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solution to a special case with unequal service times. 

Chapter Four mentioned that such a verification was attempted 
but abandoned primarily due to time constraints. It is 
quite possible that additional work in this area could be 
successful . 

The simulation algorithms should be studied to determine 
confidence limits on its estimates and to determine the 
sensitivity of the algorithm to changes in system parameters. 
This type of analysis would also serve to increase the con- 
fidence in the comparison between the simulation and 
approximation models. 

Finally, there is obviously much room for improvement 
of the analytical approximation. Although the model developed 
in Chapter Four appears to perform fairly well and has an 
intuitive rationale, there exist a number of techniques 
which could be expected to yield better results. In particu- 
lar the theory of Gaussian approximations and diffusion 
approximations discussed by Gaver and Lehoczky [Ref. 2 and 3] 
and Kleinrock [Ref. 6] should work well at least in the 
case where the system is heavily saturated. 



57 



APPENDIX A 



DETAILS OF THE DERIVATION OF THE EXPRESSION 
FOR THE EXPECTED NUMBER IN THE QUEUE (E [Q] ) 
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